Perturbative and Numerical Analysis of Tilted Cosmological Models of Bianchi type V 
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Cosmological models of Bianchi type V and I containing a perfect fluid with a linear equation of 
state plus cosmological constant are investigated using a tetrad approach where our variables are the 
Riemann tensor, the Ricci rotation coefficients and a subset of the tetrad vector components. This 
set, in the following called S, describes a spacetime when its elements are constrained by certain 
integrability conditions and due to a theorem by Cartan this set gives a complete local description 
of the spacetime. The system obtained by imposing the integrability conditions and Einstein's 
equations can be reduced to an integrable system of five coupled first order ordinary differential 
equations. The general solution is tilted and describes a fluid with expansion, shear and vorticity. 
With the help of standard bases for Bianchi V and I the full line element is found in terms of 
the elements in S. We then construct the solutions to the linearized equations around the open 
Friedmann model. The full system is also studied numerically and the perturbative solutions agree 
well with the numerical ones in the appropriate domains. We also give some examples of numerical 
solutions in the non-perturbative regime. 

PACS:numbers: 04.20.-q, 98.80.-k 



I. INTRODUCTION 

In this paper general Bianchi types V and I perfect 
fluids with linear equation of state and cosmological con- 
stant are studied. In general these spacetimes are tilted 
and in particular there are solutions with rotating mat- 
ter. It has been difficult to find exact solutions with 
both expansion and nonzero rotation of the matter flow. 
To our knowledge the only known exact homogeneous 
perfect fluid solution with rotation and expansion is the 
self-similar radiation-filled Bianchi type VIo found by 
Rosquist [1]. 

A number of rotating imperfect fluid solutions with 
heatflow are known, see e.g. [2]. Since for rotating mat- 
ter the hypersurfaces of homogeneity are tilted with re- 
spect to the restframe of matter, local space will not look 
homogeneous. Hence heatflow is expected and for some 
solutions the heatflow can be related to a temperature 
gradient [3] , but often with unrealistic coefficient of con- 
ductivity. Normally the heat conductivity is negligibly 
small and a perfect fluid approximation should work well. 

For treatments of the properties of homogeneous rotat- 
ing models in general see [4,5] and for some perturbative 
calculations see [6-8]. In [9] the qualitative behaviour 
of locally rotationally symmetric (LRS) Bianchi V solu- 
tions is analysed and in particular expressions for differ- 
ent quantities at early and late times are given. There are 
a number of works on tilted Bianchi cosmologies using a 
dynamical system approach. The irrotational subcase of 
type V was studied in [10]. Recently the late time be- 
haviour of tilted Bianchi models including type V was 
considered in [11]. The stability of non-tilted Bianchi 
models against tilt was studied in [12]. 

To find the solutions we use a method for construction 



of solutions to Einstein's equations, [13-16], based on the 
invariant classification scheme by Cartan and Karlhcdc 
[17,18]. The method is shortly described in section II. 
In section III the method is applied to Bianchi V and 
I models. First choice of frame and the set (called S) 
of quantities needed to specify the spacetime are given. 
Then the structure of the isometry group is imposed, giv- 
ing relations among the elements in S. Next the integra- 
bility conditions for the set S to describe a geometry are 
imposed together with Einstein's equations. The general 
system is reduced to an integrable system of five coupled 
first order ordinary differential equations. With the help 
of standard bases for Bianchi V and I the full line-clement 
is found up to quadratures in terms of the elements in S. 
The subclass of orthogonal solutions is easily solved, but 
all these solutions are well known. For LRS dust the 
equations can also be integrated [19,20]. 

In section IV we consider first order perturbations 
around the open Friedmann universe. The general first 
order solution depends on five constants of integration, 
the same number as for the general exact solution, and 
has nonzero expansion, rotation and shear. Finally, in 
section V the general system is studied numerically and 
the results agree well with the perturbative calculations 
in the allowed range. Examples of numerical solutions in 
the non-perturbative regime are given. 



II. CONSTRUCTION OF SOLUTIONS TO 
EINSTEIN'S EQUATIONS IN TERMS OF 
CURVATURE INVARIANTS 

A brief summary of the method is given here. For more 
details see [16,21]. 
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According to a theorem by Cartan spacetimes are lo- 
cally completely determined by a set consisting of the 
components of the Riemann tensor and a finite num- 
ber of its covariant derivatives in a frame with constant 
metric rjij [17,18]. This set will henceforth be called 

R p+1 = {Rijkl,Rijkl;m 1 ,- — Rijkl;m 1 ....m p+1 }, where p IS 

such that the components in R p+1 are functionally de- 
pendent on those in RP (as functions of both the coor- 
dinates x^ on the manifold and the parameters of the 
Lorentz group £ T ). This description can be used to clas- 
sify geometries and also to construct new solutions to 
Einstein's equations in terms of the set R p+1 . 

Assume that we have symmetries such that R p+1 
only depends on x a , a = 1,2,..., I < n =dimension 
of spacetime, (in some canonical coordinates) and rota- 
tions in the afr-planes , {";,} = l,...,m < n(n — l)/2 
(with frames adopted to the rotational symmetries). 
Here I = n — dim(orbits) and m = n(n — l)/2 — 
dim(isotropy group). The set R p+1 is completely deter- 
mined by the smaller set S = {R p q ki,-f a bi ,x a \i} where 
= Xi(x a ) = X i fl dx a /dx^ are the derivatives with 

respect to the frame vectors and Yjk are the Ricci ro- 
tation coefficients. Here the numbering is such that 
{ p q } = m+ 1, n(n — l)/2 are the complementary rota- 
tions to the { a b} ones, i.e., those that keep the set R p+1 
unchanged. 

A set R p+1 together with a constant frame metric r]ij 
describes a geometry iff certain integrability conditions, 
being equivalent to the Ricci identities (including the 
commutators for the essential coordinates) and part of 
the Bianchi identities, are satisfied. In practice it is of- 
ten easier to work in a fixed frame than to let the compo- 
nents in R p+1 depend explicitly on the parameters £ T of 
the orthogonal group. The Ricci identies then split into 
the commutators for the essential coordinates x a 

x a [li ^ lj] +x a lm7 m [ij] =0, (1) 

and the Riemann equations for rotations in the a6-planes 

R a bij = ^ a blj , a x a {i] + 27° fe b -76fei] + Zl\ k l k M (2) 

(the antisymmctrizations are only over if). When using 
the description in terms of the smaller set S equations (2) 
merely serve to define the ii a ^-components of the Rie- 
mann tensor. Note that equations (2) give 36 equations 
in the case without isotropics and the restriction to 20 
independent components come from the cyclic identity 
below or equations (1). Since not all commutors or Rie- 
mann equations are used when the spacetime has sym- 
metries some more integrability conditions arc needed. 
They are parts of the cyclic and Bianchi identities 

R\ijk] =0, t = l + l,...,n, (3) 
R p q [ir,k]=0, rj = m+l,..,n(n-l)/2, (4) 

where t — I + 1, ...,n numbers the symmetry directions 
(in a suitable numbering of the frame vectors). 



The above description in terms of R p+1 can be used 
to find new solutions to Einstein's equations. This ap- 
proach has the advantage of being coordinate invariant. 
Also, since the components in R p+l correspond to di- 
rectly measurable quantities, physical requirements arc 
easy to impose. Another advantage is that the set of 
equations (1), (2), (3) and (4) is a subset of the ones used 
in different tetrad methods like the ones in, e.g., [22] or 
the Newman-Penrose method [23]. The procedure is: 

1. Choose a set R p+1 (or equivalently S). Some of the 
elements in R p or functions of them are chosen as 
coordinates. Einstein's equations are imposed by 
restricting the Ricci tensor. 

2. Impose the integrability conditions (1), (2), (3) and 
(4), leading to a set of first order differential equa- 
tions (together with algebraic constraints) for the 
elements in R p+1 (S). 

Solving this set of equations gives R p+1 and hence a com- 
plete local description of the geometry. If the geometry 
does not have any symmetries, one can solve for all the 
1-forms from dx^ = x^lu 1 and hence get the full line- 
element d 2 S — IJijlV 1 ^ . 

When there are symmetries one only obtains part of 
the 1-forms, but one may determine the Lie-algebra of 
the isometry group, [24], and from this it is often possible 
to make an ansatz for the remaining 1-forms. Cartan's 
equations look the same in F(M), the bundle of frames 
on M, as in M 

dJ =u> j AwV (5) 
dw i j = -w i k *w k j + ±I? jld u k *u l (6) 

but now d = d x + d%. The connection forms uj 1 - will 
now be linearly independent of the 1-forms a/, and can 
be written as 

lo) = 7 > fc + t)= 7 > fc + r) T d e (7) 

where are the generators of the orthogonal group. 
With the same notation for indices as above Cartan's 
equations split into those corresponding to essential coor- 
dinates and rotations (x a and r%) and those correspond- 
ing to symmetry coordinates and rotations {x t and r p q ) 

du* = ujj Aw'j, duj\ = -uj p k ALu k q + \R p qkl Lu k f^J . (8) 

The 1-forms that arc missing can be found from these 
equations. Since spacetime locally is completely deter- 
mined by R p+1 , it is sufficient to find one linearly in- 
dependent solution to equations (8). All other solutions 
will be locally equivalent to this one up to coordinate 
transformations. When projecting onto the orbits, given 
by dx a — and r a b — 0, equations (8) reduce to the 
Lie-algebra of the isometry group (since the group acts 
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simply transitive on the orbits in F(M) [25]) and the 
structure constants will be given by the elements in R p+1 
[24]. If one wants to find the full line-element, the pro- 
cedure can therefore be continued in the following way 

3. Solve for the u a (in a suitable numbering) and the 
Lu a b from R p+1 . 

4. Determine the Lie algebra of the isometry group 
and some standard 1-forms satisfying it. 

5. From these make an ansatz Co 1 for the 1-forms over 
the full spacetime. 

6. Calculate the corresponding R p+1 and compare 
with R p+l , giving differential or algebraic equa- 
tions for the coefficients in the metric ansatz. These 
equations are of course equivalent to (8). 



III. BIANCHI V AND I 

In this section we consider homogeneous cosmological 
models of Bianchi type V and I, i.e. those character- 
ized by the symmetric matrix in the Ellis-MacCallum 
scheme [22] being zero. We assume that matter can be 
described as a perfect fluid. First the preliminaries, like 
choice of frame and the elements in the set S (R p+1 ) 
are given. Einstein's equations with a cosmological con- 
stant are used. Then it is illustrated how one imposes 
the isometry group. After this we give the integrability 
conditions and reduce them to an integrable system of 
five first order ordinary differential equations. 



A. Preliminaries 

As energy- momentum tensor we take that of a perfect 
fluid 



(p + p)uiUj -prjij 



with linear equation of state p — (7 — I) p. Here p is the 
restframe density, p is the isotropic pressure and u % the 
4-velocity of matter. Since homogeneity is assumed the 
elements in the set S will depend on only one timelike 
coordinate, that we choose as the density p. Sometimes, 
especially in problems with more than one independent 
coordinate, it can be advantageous to specify the coordi- 
nates at a later stage to simplify the equations, see [21]. 

A Lorentz frame u> 1 will be used. We choose a comoving 
frame, i.e., the 4-velocity is given by u — Sfu; 1 = ui°. In 
general the normals of the hypersurfaces, dp = p\iUj l , 
will be tilted relative to the 4-velocity. The 1-direction 
is chosen to be in the direction of the spatial part of 
the density gradient, i.e., p\ 2 = p\3 = 0. From these 



equations we see that, once p\ and p\i are determined, 

one can solve for uj° as 

o dp p\i j d 

uj' = lu or X Q = p\ —. 

P\o P\o dp 

This choice of frame means that we deviate from the 
usual approach of adopting the frame to the hypersur- 
faces of homogeneity. 

Since there is only one essential coordinate, p, this is 
the only 1-form that we will be able to solve for from 
R p+1 . The frame is finally fixed by requiring that the 
vorticity (rotation) vector of the fluid is in the 12-planc, 
i.e., fi 3 = \e z% i k u)ijUk = 0, corresponding to that U12 = 
(see below for the definition of the vorticity tensor). 

The general model in this class will not have any 
isotropies and one can from the set consisting of p\ , p\i 
and all -f l jk construct the full set R p+1 . (The set could 
in principle be even more reduced since the equations (1) 
give relations among the quantities.) Since we want to 
impose Einstein's equations 

Gij — ^ij ~t~ AT^ij 

for a perfect fluid and are using a comoving frame the 
Einstein tensor should be given by 

Goo = p + A, Gn = G 22 = G 33 = p - A, 

where A is the cosmological constant. The cyclic iden- 
tity must also be imposed. The nonzero elements of the 
Ricmann tensor are then given by 



R 



0101 



Gi 



1 1 

2 P ~6 P 



1 



A, R 



•0102 



R 



1323 



G 2 , 



#0103 — — #1223 — G3, #0112 — #0323 — G4, 
#0113 = —#0223 = G5, #)123 = Cq, 

#0202 = G 7 - -p — -p + - A, i? 203 = #1213 = G 8 , 

2 6 3 



#0212 

#0303 



"#0313 — Gg, #)213 — G 



1 



1 



-C 1 -C 7 --p--p- 



10, 



:A, R, 



0312 



G10 — G 



6, 



#1212 — Gi + G7 — -p 



-A, R 



1313 



#: 



2323 



-Gr 



1 X A 

3^-3 A 



(9) 



where Ci are the ten independent components of the 
Weyl tensor. 

Some of the rotation coefficients are expressible in 
terms of the kinematic quantities shear , vorticity uiij , 
expansion 9 and acceleration a; t 



-LUij — (Tij + —hijO — CLiUj. 



where a tj = hfh/ (u^.q + \h kl 6), uj i:j = tifh/ 
8 = u 1 ^ and = u^jV? and hij — UiUj — r/ij is the 
projection operator onto the space perpendicular to the 
4-velocity. In the next subsection we impose the require- 
ment that the isometry group is of Bianchi type V or I. 
This will give six restrictions on the rotation coefficients. 
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B. Symmetry group 



C. Integrability conditions 



The Lie algebra of the isometry group can be deter- 
mined by projection of Cartan's equations onto the orbits 
[24] . Since we do not have any isotropies in this case, the 
orbits will be the hypersurfaces of homogeneity dp = 
in M. From dp = we get 



dp = p| w°| + p\ a u: a \ = or 



■to 



where a — 1,2,3 are the spatial indices and a vertical bar 
| indicates projection onto dp = 0. From the requirement 
of no isotropy one has that r*-| = holds on the orbits 
and (7) then gives 

w j\ = 7 jk^ I • 

Hence the orbits are spanned by {w |,w |,w |}. By pro- 
jecting the first pair of Cartan's equations (5) (the second 
pair will not give anything new in this case) one obtains 

duj a \ = a;- 5 1 A = ^C a ps u p \ A u 5 \ 
where the structure constants are given by 

2 C es- 1 m + —^m-—^m- ( 10 ) 

The structure constants are hence given in terms of ele- 
ments in R p+1 , and since they are functions only of p the 
Cpg are constants on the orbits. 

From the Ellis-MacCallum scheme for the Bianchi 
classes, [25], we can decompose the structure constants 
into one symmetric 3x3 matrix N a @ and a 3-vector A a 
according to 

\c% s e^ = N a ~< + e a ^A giving 



2~ 06 



(11) 



Bianchi classes V and I are characterized by N aS = 0. 
This will give six relations among the 7* fe . Note that 
we cannot simply use the standard form of the struc- 
ture constants in (10) and obtain nine relations for the 
Ricci coefficients this way since our choice of frame (giv- 
ing other relations among the elements in S) might be 
inconsistent with that giving the structure constants in 
standard form. The nonzero Ricci rotation coefficients 
are then given by (a, (3 = 1, 2, 3 and uj\2 = 0) 



70a/3 = -<Jaf3 ~ U a f3 (a ^ /?), 



7oaO = — a Q , 

70aa = \6 - &aa, 7123 = 7132 = — ^<T23, 



7133 = 7122 + (fH + 2(T 2 2) , 7231 = (7230 - ^23) j 

7232 = 7131 + (fl3 + ^13 - 7130 ) , 7130 , 7131, 7230, 

7233 = -7121 + (— fl2 + 7120) , 7120, 7121, 7122- 



The commutator equations (1) give 

^i-^i + P|*(^-7*«) = (13) 
and the Riemann equations (2) 

Rl 3 ki = n ]%P P\u\ + 27 im [i7 jro fc] + 27^7"^] (14) 

(antisymmetrization only over fcZ). The cyclic identity 
is already imposed due to the choice (9) of the Riemann 
tensor and the Bianchi identities need not be imposed due 
to the lack of isotropies. (However it is in practice often 
convenient to use the twice contracted Bianchi identities 
since they give simple relations, but they may of course 
be derived from the other equations). Hence the system 
(13, 14) is the complete system. It is a set of first order 
ordinary differential equations, where the independent 
variable is p, and algebraic constraints. Some are easily 
solved and give: 

0-2 = CI3 = 7121 = ^23 = 0, 7120 = 012, 7130 = 013 +^13, 

7230 = -o-23 (or 7i3i = 0), p\ = - jp9, ai = (1 - j)v9, 



t*>13 = 2^31 

where we have introduced the tilt 
v = P\i 
P\o ' 



(15) 
(16) 



The case 7131 = will be treated separately in section 
IIIF. We see that the vorticity given by 

W13 = ^7131 or n 2 = -7^7131 (17) 

is perpendicular to the acceleration, a result already 
found in [4]. The ten components, d, of the Weyl tensor 
are also given by (14). 

The system then reduces to nine first order differential 
equations 



fi(p)=Fi(fM,p), 
where / = df /dp, for the functions fi(p) 

fi E {v,0, an, 022, 012,(713, C723, 7131, B} , 

where we use the variable 

B = 7122 - \v0 + v(T 2 2 



(18) 



(19) 



instead of 7122, and four algebraic constraints 

G a (f j (p),p)=0. (20) 



(12) 



It turns out that the differentiated constraints, when us- 
ing (18), all are satisfied, i.e., 

dG a 

Of/ 3 ' dp ~ dfj' 3 • dp 
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Hence the system can be reduced to five differential equa- 
tions and 4 constraints [15]. The constraints are given by 

7i3i [~2v 2 (3(7n + (3 7 - 4)0) - ivB- 
18 (on + o 22 ) (l - v 2 )] - 18Bcri3 = 0, 

Bu 12 - 7131C23 (l - v 2 ) = 0, 

97131 (7i3i« + 2o-i3 ) (1 - v 2 ) + 6v (2B 2 + 7 p) 
-2v 2 B [2(3 7 - 4)6> - 3ou] - 18Bau = 

97131 [7131 (12 - 5v 2 ) - 8vai 3 ] + 48vB9 + 



36 [a 2 u 



O11O22 + + of 2 + O13 + 023 + A + 3B 2 



+p - v 2 (o n 22 + 022 + + 402 [( 6 7 - 5)« 2 - 3] 



+ 12oii [<ovB + (3 7 - 2) v 2 9] = 



(21) 



that are easily solved for 013, 023 and on in which the 
three first constraints are linear. Finally the last con- 
straint gives a second order polynomial in o 2 2- The dif- 
ferential equations are 



v = v 



± _\1_ 

IPO \ 3 7 / P 



7131 = 7131 



012 = 012 



B = B 



(011 + 022) 

ipQ 

(011 — 022) 

ipQ 

011 1 

7p0 3jp 



+ 



1 

37/5 

J_\ _ 2^7131023 

IP J ipo 

27131013 
IP 



(9 2 J = [37i 2 3 iw 2 + 12(7 - l)v6 (B + van) - 20 2 + 
2(7 - 1)(6 7 - b)v 2 9 2 - 3 ((3 7 - 2)p - 2A) - 

■4s)] 



O11O22 + 012 



12 (0 

/[3 7 p((7-l)« 2 -l)] 



'13 



'22 



(22) 



If one instead wants to use proper time for a comoving 
observer as independent coordinate the equation 



dp 



P\o 



-7,30. 



(23) 



obtained by putting dx — dy — dz = in the linc-elemont 
(see equation (26) below), should be added to the system 
(22). The derivatives in (22) are then expressed as 

I = djj_ = dh dr = dfj 1 
* dp dr dp dr jp9 

The system has a unique solution for given initial con- 
ditions provided for the reduced system (22) satisfy a 
Lipschitz condition (for example if Fj are C 1 in a compact 
convex domain). The general solution hence depends on 
five constants of integration and in general its matter flow 
has both expansion, shear and vorticity. 



One could of course choose to solve for four other 
quantities than 013,023,011 and 022 from (21). When 
doing a first order perturbative calculation around an 
isotropic universe, as in section IV, both on and 022 
cannot be solved for from (21) since the last constraint 
will be quadratic in o 2 2 or on. Hence we solve for 9 from 
the last constraint and use the differential equation for 
022 

022 = °- + [(7 - 2)p - 4 (o 2 2 - (v 2 - 1) o 2 3 ) - 2A+ 



2 022 



2(7n + (67-5): 



AB 2 



3/ V V 3 , 
+2vB (2(022 - on) - l9) + 47131 (uoi3 - 7131) ] 



/ [2 1P 9{v 2 1)] 



(24) 



instead of the one for 9. 

This system is not in a suitable form for a dynami- 
cal system analysis. The system should then be made 
autonomous and compact dimensionless variables should 
be introduced. In [10] irrotational tilted Bianchi type 

V cosmologies were studied with this method. The field 
equations are derived in terms of expansion-normalized 
variables making the state space compact. The existence 
of a monotonic function shows that the dynamics to a 
large extent is determined by the invariant subset of lo- 
cally rotationally symmetric models. A complete analysis 
of the orbits with non-extreme tilt was obtained. In [11] 
tilted Bianchi models of solvable type, including Bianchi 
type V, were considered with emphasis on the late-time 
behaviour. The equilibrium points were given in [10], but 
here the stability analysis was performed in the full state 
space. It was found that for 2/3 < 7 < 2 Bianchi type 

V models approach the Milne universe in the asymptotic 
future. To complete the analysis a study of the behaviour 
for early times (high densities) would be of interest. 



D. The metric 

From equations (10) for the structure constants we find 
the Lie-algebra of the isometry group to be 

(Lo 1 \ = 7131W 1 ! A oj 3 \ 
duj 2 \ = Bu 1 \ Aw 2 | + 7i3iw 2 | A iv 3 \ 
duj 3 \ = Buj 1 ] A u) 3 \ 



(25) 



where B is given by (19). As expected it is not in the 
standard form for Bianchi V, unless B = 0. Guided by 
the above algebra and standard 1-forms for Bianchi V, 
we make the following simplified, but sufficient, ansatz 
for a basis of 1-forms for the full spacetime 



5 



OJ 



Co 2 = g\e x dy + g 2 e x dz + g 3 dx 

1 B 

dx + hi e x dy 

7131 7131 



dp -1 -1 , -x , 

- — n~ vu} . w = /i e «y 



£ 3 = 



(26) 



where /i, g\, g 2l gz and /ii are functions of p to be deter- 
mined. (The case 7131 = is treated separately in section 
IIIF.) From (26) we calculate the set S — {p~u, Jijk}- A 
comparison with the set S gives 

/1 = hi = Cipy^v , 31 = 7i3iP 1_ ^^i 

.92 = C 2 7i3i/? 1 ~"v , ff 3 = 7i3i/> 1 ~~^2 (27) 

where C\ and C2 are (nonzero) constants of integration 
that can be absorbed in the definitions of the coordinates 
y and z and the integrals Ii and I 2 are given by 



/1 = 



2Ci 



^23/3 



2 /• 

7 J vpOjl 



v 2 p0-f( 3 i 
1 

-dp 



[0-127131 + 023-B] dp 



(28) 



respectively. From the 1-forms the metric is given by 



The full metric is hence given in terms 



of quadratures once the set S has been constructed. The 
solution depends on some arbitrary constants of integra- 
tion and an even more general ansatz than (26) could 
have been made, but all these metrics give the same set 
S and hence they are locally equivalent. 

Later on we will consider perturbations around the 
open Fricdmann model, and will be interested in the limit 
v — > as e — > 0. As seen some of the components in the 
metric will then diverge. This coordinate singularity can 
be avoided by changing coordinates according to 



y = ty 



z = 



E. Orthogonal solutions 

For the orthogonal case v = 0, when also the vorticity 
is zero, it is better to use a frame where the shear tensor 
is diagonal, <Ti 2 = 013 = 023 = 0. The obtained system is 
easily solved and all solutions are well known. Essentially 
only two types of solutions appear. These are the Bianchi 
V solutions 



on 



-kipi , 7131 = 7232 = k 2 p 3 -i , 

6 = ±V3\J 3fc| p^ + k\p$ +p + A (29) 
and the Bianchi class I solutions 
on = hp-y , cr 22 = k 2 p-< 
9 = ±VW (kf + k 2 + kik 2 )p^ +p + A . (30) 



All orthogonal solutions can be found from these inter- 
changing the 1-, 2- and 3-directions. 
The corresponding metrics are 



dp 
7P#' 



cie 



and 



c 2 e 



dp 

_ r T022 

c 2 e 3 ~<i> e 



Ldp e- X dz, uj 3 



—p 3 -<dx 
k 2 



lo 1 = c x e $ y °p*> dp dx , 



dp dy, 



c 3 e 



respectively. The integrals can be performed for specific 
values of 7. 



F. Irrotational tilted solutions 

In [10] this class for non-extreme tilt (v < 1) was stud- 
ied using a dynamical systems approach. It was found 
that that the dynamics to a large extent is determined 
by the invariant subset of locally rotationally symetric 
(LRS) models. This subset was studied in detail in [9] 
for different values of 7. There are solutions for which the 
tilt always is less than one and hence the hypersurfaces 
of homogeneity remain spacelike when going backwards 
in time. However, it may approach one for large times for 
certain values of 7 including the radiation case 7 = 4/3. 
There are also solutions for which the tilt gets larger than 
one for small times and hence the hypersurfaces change 
from being spacelike to being timelike. For this class sin- 
gularities may occur for finite densities. For example, in 
the 7 = 4/3 case some of the kinematic quantities and 
the Weyl tensor diverge for v = V3. 

The irrotational solutions are given by 7131 = 0, which 
implies that the vorticity vanishes and the frame cannot 
be fixed by demanding O 3 = 0. Instead the frame is fixed 
by putting 023 = 0. Some care should be taken in using 
equations (21) and (22) directly since an extra constraint 
is introduced and they were derived assuming 7131 7^ 0. 
Yet the only nontrivial cases are those obtained from this 
system with 7131 = 023 = and B / 0. From the two 
first constraints eri 2 = 013 = is obtained, an and 022 
can be solved for from the two others. The system of 
differential equations is reduced to a system for v, B and 
9. This can be reduced to a system of two differential 
equations since the equations for v and B in this case 
can be combined to 



i = § + (I-i)i 

v B V7 J P 



with solution 



B = Civp- 



(31) 



6 



where C\ is a constant of integration. The metric is in 
this case given by: 



giving 



VLJ 1 , LU 1 



7p(9 

i-) If cr 22 j 

w = k 2 p~^e J ~ a p e~ x dy, 

uj* = k 3 p 3 ->e ■' ^» p e ^dz. 



L MS solutions 

LRS solutions are obtained by putting cru = — 2(722 
in the above equations (the LRS symmetry lies in the 
23-planc). If one solves the constraints for CT22 and 9 the 
system is reduced to one differential equation for v. 

For 7=1 and A = the equations can be completely 
integrated [19,20]. The solution in our variables is given 

by 



(T11 — ~ 2(722 = i 



K 2 (3CK + 2) 



v = T 



3(D- «V1 + Ck) ' 
3£W1 + Ck -k 2 (2+ |Ck) 

7122 = -K 



(32) 



where C and D are constants of integration and k and p 
are related as 



3CDk 3 



D - ksJI + Ck 
The basis 1-forms are 
dn 



(33) 



2 7 3 7 

w = ay , w = az . 



=F ax , w = ax , 



\ Vma *\- 1 + |Ck 2 

that is less then one and hence the hypersurfaces of ho- 
mogeneity remain spacelike for all times. 

A perturbative calculation (see IV for the perturbative 

method) with a small 7131 = 67131 around the exact 
solution gives 



7i3i ^ F 



D - kVT+Ck 



(35) 



(34) 



where F is a constant of integration. Hence the per- 
turbation grows as the other quantities for the first case 
(C,D > 0) as k — > Ki. However, for the second case 
(D < 0) it goes to zero when k — > 00. 



G. Solutions with a timelike homothetic motion 

A spacetime has a homothetic Killing vector £ if 

is satisfied for some constant c. This implies that quan- 
tities of the same dimension scale in the same way in the 
direction of £. For a timelike homothetic Killing vector 
it is hence easy to show that for perfect fluids all compo- 
nents of the Riemann tensor, Rijki , are proportional to 
p, the Ricci rotation coefficients to p 1 / 2 and p^ to p 3 / 2 , 
see e.g. [26]. The field equations are then reduced to 
algebraic equations. 

Unfortunately this limitation is quite restrictive and 
the only solution of this type is the flat Friedmann uni- 
verse (of Bianchi type I) . 



IV. PERTURBATIVE SOLUTIONS 



From (33) one finds essentially two types of solutions. 
If both C and D are positive the density rises from zero 
to infinity when k goes from zero to Ki, where Ki is given 

by 

D - /ci-v/i + CkI =0. 

The tilt v then also increases from zero to infinity, and 
hence the hypersurfaces of homogeneity change from be- 
ing spacelike to being timclikc. 

A positive density is also obtained for C > and D < 
0, and in this case the density goes from zero to infinity 
when k goes from zero to infinity, but now the tilt varies 
from zero to zero in this interval. The maximum value is 
obtained for k 2 



CkI + 2D^/1 + Ck 2 = 0. 



In this section we consider perturbative solutions to the 
system (21) and (22). In general perturbative solutions 
to systems of non-linear equations with constraints need 
not correspond to exact solutions, but in the present case 
the system can be reduced to (22), that under reasonable 
assumptions is known to be intcgrable with five constants 
of integration. Hence a perturbative solution, obtained 
by solving the Taylor expanded equations, should agree 
with the Taylor expansion of some exact solution. This 
is also favoured by comparing the perturbative solutions 
with a numerical solving of the full system, see section 
V. 

We here first briefly recall some basic results of per- 
turbation theory. Calling the four functions solved for 
from the algebraic constraints (20) g a , a,b,.. — 1,...,4, 
the reduced system can be written as 
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fi = Fi(fj,g a ,p), = 1,...,5, 

= G a (f j ,g b ,p). (36) 

The functions fi and g a are expanded around the solution 
and g^ of (36) in the small parameter e as 



Ai^-./f ) = £ if ) + e 2 ./f ) + ... 
Ag a ^g a -g^=eg^+e*gV + ... 

giving, with H p equal to Fi or G a , 
^ = Wf, 5 i°),p) + ^A/ J + ^A, a . 

d 2 H p » j» » 1 <9 2 iJ„ a j* a j? 
dfjdg a Zdfjdfk 



(37) 



i <9 2 #,, 



^A 5o A ffb + ... = J ff p (.f (0) , ff (°),p) + 



2 



df: 



9a 



1 <9 2 ff 



2dfjdf k J > 



2 dg a dg b 



dfjdg a 



E-fW£)) + ... (38) 



(0) 



where the partial derivatives are evaluated at fi = f % 
and g a — g^ . Identifying equal powers of e in (36), 
using (37) and (38) with H p = G a , one can solve for g^ 



in terms of f- n ^ and lower order quantities as 



(39) 



where (G x ) a is the inverse of (assuming that it is 

invertable) and (? a ™ -1 ^ only depends on f- m ^ and g^ 
up to order n — 1, i.e. m < n — 1. Substitution of (39), 
(37) and (38) with H p = Fi into (36) and identification 
of equal powers of e now gives 

,(„) / dF, dFi , nt dG b \ ( „) ( „_i) 



where J 7 - 11 only depends on f- m> and g a " l ' up to order 



n — 1 and for n = 1 reduces to J 7 , 



(n-l) 



= 0. From 



this equation we recall the result that the homogeneous 
parts of the differential equations are the same to each 
order in e, and hence the integration constants add up 
as ki — ekj 1 ^ + e 2 fc| 2 ' + so that the correct number of 
constants is maintained to any order. 

A. First order perturbations 

We here consider perturbations to first order in the 
small parameter e around the open Fricdmann universe. 



The nonzero elements in the set S for the Friedmann 
universe are given by 



g(0) 



±V3W3(fc* 0) ) p^+p + A, 



(o) _ (o) 

7l33 — 7l22> 



(0) 
7131 



where fcf^ is a constant of integration. In the following 

we choose 7^2 = 0, so that 7^ = 7232 = kf^p^. The 
freedom in one of the Ricci rotation coefficients is due to 
that only the 7oij appear in S for the Friedmann universe 
due to the isotropy. Note, however, that the resulting 
perturbed solutions are depending on this choice. For 
example, if we instead had chosen 7^ = 0, the first 
order perturbations would have had zero vorticity. 

Instead of solving for 022 from the last of equations (21) 
we solve this equation for 9, and hence the differential 
equation for 9 in (22) is replaced by the corresponding 
for cr 2 2, (24). 

To first order we write the elements in S as 

(1) (o) , (1) 

Expanding the system (21) and (22) (with the last equa- 
tion replaced by (24)) to first order then gives 



(o) 

7232 



± 



(k^y 



(41) 



)3 7 





= ek 3 9p : ^ 1 




= -022 = efc 




7&3 I 


Cl3 


3fei 




k^ 3 5 _ 


Wl3 


= ' 2 ^ 



7122 



£7l22 

k 3 J_ 
-TP 31 



where 
2 



7122 



7131 = fcl/0 3T : 
1 

k 5 p 3 ~i - 



7 



(7-1)^ 



P" 



e 



dp 



9 = ±V3J3kfp 



where fci = k\ 



(0) 



ek[ 1} , k 2 



•A, 



k£\ k 3 



(42) 



and fc 5 = k\ are the five constants of integration. 

The integral in 7122 can be evaluated for certain values 
of 7 if A = 0. For 7=1 (dust) one gets 



(1) , i , foo 

7l22 = hp 3 ± — 



(pi -6k!) , 



for 7 = I (radiation) 
J122 = hp* - 



hp' 4 



-9p * ± fci In 
9 F 



'pi± ±6p~i -V3h' 
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and for 7 = 2 (stiff matter), with fc 5 = k§ ± k 3 /2, 



7i22 = hp'T 



4 In 



± 



5 e \ 9\/3fc?p-5 



V^fci 3fei / V3p^ ± 6» 



The integral can also be evaluated with nonzero A for 
7 = 1 if ki = 0, corresponding to that the background 
metric is of Bianchi type I (but the perturbed one is of 
type V), giving 



7& 



; i &3 I / — T 



B. The metric to first order 

The form of the metric (26-28) is not suitable for a 
perturbative calculation since some of the coefficients di- 
verge when e — > 0. This can be avoided by introducing 
new coordinates according to 



y = ey, z = - . 

e 



(43) 



When expanding the 1-forms to first order in e the tilt v 
to second order, given in IV C, will be needed. To first 
order in e one then obtains the following 1-forms (with 
Ci = C 2 = 1) 

dp a 1 

w u = — - ep-i L e x dy 

7P<9 

= -^-p~^e~ x (1 - eAi) 
w = -e— -Aip 3 7e 

hk 3 p~^e~ x (1 7 eAi) <iz 

3 1 , 

or = — p .^ dx 7 

fci 

^-(Ip-^-ip*- 1 ^ (44) 
where and 7^2 are given by (42) and A\ by 

A 1 (p) = -k 2 J P —^dp . (45) 

For dust (7 = 1) and A = A\ becomes 
2k 2 9 



. 2k 2 V ( 2 _i\ 

M = — — (l-6fc?p 



and with A^O and k\ = 



2fc 2 6> 



Radiation (7 = 4/3) with A = gives 

Ai = -k 2 6p- 1/2 . 

The relation to proper time for a comoving observer is 
obtained by putting dx — dy = dz = in the line-element 
and is given by dr = ^. For example, for an expand- 
ing universe with 7 = 1 we have 9 = V%\J Zk 2 pi 7 p. 
Integration gives (with integration constant such that 
r(p = 00) = 0) 



6p 



1 



In 



3kip3 



3k 2 6kf \e-3k lP i J ' 

Note that deviation from zeroth order quantities in the 
expression for proper time will appear first to second or- 
der. 



C. Second order perturbations 

From (39) we see that the second order perturbations 
can be obtained up to quadratures. A full second order 
calculation will be done in a forthcoming paper. Here we 
focus on the tilt, v, whose n:th order equations decouple 
from other n:th order quantities, and also since it will 
be needed to second order to obtain the metric to first 
order. 

To second order one obtains 




where now fc 3 = k 3 ' 7 ek 3 ' . For A = and 7 = 1 we 
have 



v = ek 3 p? ( l-^f 2 ^i ^2 



1 - 6kfp-*) 
ek 3 pi (l T 2£|^,*+3fc?(,*-6fc?)) 
and for A 7^ 0, 7 = 1 and k\ = 

v = ek 3 p^ (l 7 + A ) • 

For A = and 7 = 4/3 the tilt is given by 
v = efc 3 ^1 — ek 2 9p~i^J = 

eh (l T ek 2 y/3yj + 3kfj . 

These results are in accordance with the exact equation 
(22) for v that shows that the behaviour is determined 
by the sign of an/9 (an — efc 2 p 1 / 7 ). 
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V. NUMERICAL SOLUTIONS 

In this section we solve the system (21-22) numerically 
using the Runge-Kutta method with a truncation error 
of order h 4 in the step length h. A comparison between 
the code and the exact solutions (29) gives agreement to 
high accuracy. 

A. Comparison with perturbative calculations 

To check the perturbative calculation in section IV (or 
conversely the numerical method) we have solved the sys- 
tem numerically for small deviations (at start) from the 
Fricdmann models and compared with the perturbative 
solutions. The constants of integration are chosen so that 
the numerical and perturbative solutions agree for the 
starting value of the independent variable p. In the two 
following figures first and second order perturbations to- 
gether with the numerical solutions for v with 7 = 1 (fig. 
1) and 7 = 4/3 (fig. 2) are depicted. The agreement is 
similar for the other quantities in S and other values of 
7- 

0.075 r 
0.07 - 



0.065 - 




P 



FIG. 1. Comparison between numerical calculation (solid 
curve) and perturbations to first (dashed) and second order 
(dotted) for v in the case of dust and A = 0. Initial values: 
7i3i = 1, v = 0.03, B = 0.03, ai2 = 0.03, cr 22 = 0.03. 



0.032 r 



0.0315 - 



0.031 - 



0.0305 



0.03 




0.0295 - 

0.029 I i i 1 1 1 1 1 1 1 ' 

1 23456789 10 11 

P 

FIG. 2. Comparison between numerical calculation (solid) 
and perturbations to first (dashed) and second order (dot- 
ted) for v in the case of radiation and A = 0. Initial values: 
7131 = 1, v = 0.03, B = 0.03, (Ji2 = 0.03, cr 22 = 0.03. 



B. Asymptotic behaviour 

In [11] a dynamical system analysis was used to find 
that for 2/3 < 7 < 2 Bianchi type V models approach 
the Milne universe in the asymptotic future (for low den- 
sities). Figure 3 shows the components of the shear nor- 
malized with expansion (aij/9) for 7=1, and they all 
vanish in the limit p — ► 0. The tilt v also approaches zero 
as shown in figure 4. 

0.15r 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

P 

FIG. 3. crij/8 — -> as p — » for dust and A = 0. Initial 
values: p = 0.8, v = -0.09, 7131 = 0.08, B = 0.1, a 12 = 0.11, 
CT22 = 0.12. From top to bottom (at p = 0.8) un/9, an/0, 
0-12/6, U22/6 and (T23/6 
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- 
-0.01 - 




-0.09 1 ' 1 1 1 1 1 1 ~— - — 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

P 

FIG. 4. v — ► as p — > for dust and A = 0. Initial values 
as in figure 3 

Figure 5 show the corresponding behaviour with a non- 
zero cosmological constant. 

0.1 r 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 



FIG. 5. Oij/6 — *■ as p — > for dust and A = 1. Initial 
values and ordering same as in figure 3 

As found for the LRS case in [9] anisotropics may re- 
main for late times with 7^1. This is not in conflict with 
the result of [11]. An anisotropy remains in the matter 
fields, but since matter is getting infinitely diluted space- 
time still approaches isotropy. In figure 6 W13 /9 and v are 
plotted for 7 = 4/3. As seen they do not approach zero 
for small p. 



0.02 r 



- 



-0.02 - 



-0.04 - 



-0.06 - 



-0.08 - 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

P 

FIG. 6. 0113/6 (upper curve) and v as p — > for radiation 
and A = 0. Initial values: p = 0.8, v = -0.09, 7131 = 0.08, 
B = 0.1, (T12 = 0.11, 022 = 0.12. 

For the LRS case two different behaviours close to the 
initial singularity was found in [9]. One where the tilt v 
goes to zero and one where it passes the extreme value of 
one and grows towards infinity. In the second case, how- 
ever, often singularities in some of the other quantities 
occur while v and the density still are finite. The nu- 
merical runs seem to show that generically the tilt grows 
towards one, even if it initially may be decreasing for a 
long density interval. In figures 7 and 8 v, an and 022 
for 7 = 1 and 7 = 4/3 are shown respectively. For these 
particular cases v was decreasing for as far the numeri- 
cal calculations could be carried out (about 20 times as 
long as shown in the picture). Note that asymptotically 
022 = — oil and hence these spacetimes do not approach 
LRS. 

0.6 r 



0.4 



0.2 







-0.4 




-0.6 - 



50 100 150 200 250 300 350 

P 

FIG. 7. v, an 1 and 022/ for dust and A = when 
p — » 00. Initial values: p = 1, v = 0.2, 7131 = 1, B = —0.2, 
<7i2 = 0, G22 = 0.7. From top to bottom an/6, v and 022/0. 
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6- 



0.4 



0.2 







-0.4 




-0.6 - 



50 100 150 200 250 300 350 

P 

FIG. 8. v, an/ 9 and 022 /0 for radiation and A = when 
p — > 00. Initial values: p = 1, v = 0.2, 7131 = 1, B = —0.2, 
C12 = 0, (J22 = 0.7. From top to bottom an/8, v and 022/6- 

Small changes of the initial values are sufficient to make 
v eventually turn and start growing as shown in figure 9. 



where B is given by (31). For the negative root the tilt 
may grow larger than one as seen in figure 10: 




"0 50 100 150 200 250 

P 



FIG. 10. v as a function of p for radiation and A = with 
the negative root in (47). Initial values: p — 0.8, v = 0.3. 

a = 1. 




FIG. 9. v for dust and A = when p — > 00. Initial values: 
p=l,v = 0.2, 7l3 i = 1, B = -0.2, o-i 2 = 0.01, a 22 = 0.7. 



Due to an apparent singularity in the equations for 
v = 1 we have not been able to follow the evolution 
beyond this value. For the LRS case it is possible to 
rewrite the only remaining differential equation to avoid 
this problem, but in the general case the expressions be- 
come quite complex. The equation for 7 = 4/3 in the 
LRS case is given by 

v _ (3B 2 + 2p) [2v (9B 2 + 2v 2 p) ± 



v 2 2p (3£ 2 (9 - v 2 ) + Av 2 p) {9B 2 + (v 2 + 3)p) 



± (v 2 - 3) v^l-B 4 + 9B 2 v 2 p + 27B 2 p + 4v 2 p 2 



(47) 



For this case, however, the kinematic quantities as well 
as the Weyl tensor diverge for v = y/3, as also found in 
[9]. 



VI. CONCLUSIONS 

In this paper it was shown that the general Bianchi V 
cosmology with linear equation of state and cosmological 
constant can be reduced to an integrable system of five 
ordinary first order differential equations for quantities 
that give a complete local description of the geometry. 
The full line-element was found in terms of quadratures 
of these quantities. In general the solutions have expan- 
sion, shear and vorticity. The system was cast in a form 
suitable for perturbative calculations and the first order 
perturbations around the open Friedmann model with 
vorticity, being approximations to exact solutions, were 
constructed. Perturbative calculations to higher orders 
would be straightforward up to quadratures. 

A numerical study was done and the results agree well 
with the perturbative ones in the appropriate domains. 
For large times (small densities) the results agree well 
with previous works [9-11]. Numerically we found that 
the tilt probably falls off towards zero when the density 
grows for special initial values, but generically it seems 
as if the tilt eventually always grows unlimited for large 
densities. In [10,11] the late time behaviour of (among 
others) Bianchi V solutions were studied using dynamical 
systems analysis. It would be of great interest to extend 
this analysis to early times. 
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